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(-H The conventional wisdom for the formation of the first hard binary in core collapse is that three-body interactions of single stars 

rS form many soft binaries, most of which are quickly destroyed, but eventually one of them survives. We report on direct N-body 

I simulations to test these ideas, for the first time. We find that both assumptions are often incorrect: 1) quite a few three-body 

^ interactions produce a hard binary from scratch; 2) and in many cases there are more than three bodies directly and simultaneously 

-*-^ involved in the production of the first binary. The main reason for the discrepancies is that the core of a star cluster, at the first 

^ deep collapse, contains typically only five or so stars. Therefore, the homogeneous background assumption, which still would be 

I I Reasonable for, say, 25 stars, utterly breaks down. There have been some speculations in this direction, but we demonstrate this 

result here explicitly, for the first time. 
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J.. Introduction 

. Nowadays, it is recognized that the dynamical evolution 
6f a star cluster is mostly understood as follows (Heggie 
& Hut, 2003). Two-body relaxation drives core collapse. 
The core collapse leads to high stellar density region at the 
cluster center. Owing to the high stellar density, binaries 
hie formed through three-body encounters (Aarseth, 1971; 
ileggie, 1975; Hut, 1985; Goodman & Hut, 1993). These 
binaries play important roles for post-collapse evolution, 
triggering re-expansion of the core, as well as gravothermal 
oscillations (Bettwieser & Sugimoto, 1983; Makino, 1996). 
However, it has not been confirmed whether the bina- 
ries are actually formed through three-body encounters, al- 
though it has been well investigated how they evolve subse- 
quently, after formation(Spurzem & Aarseth, 1996). Indi- 
rect evidence has been gathered indicating that more than 
three stars might be involved (McMillan, 1988), but no di- 
rect observations have been published of the core conditions 
that lead to the formation of the first binary, in a simulation 
in which a cluster of single stars undergoes core collapse. 
In contrast, later hardening of binaries has been described 
in considerable detail by McMillan et al. (1990) (MHM90). 
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In this paper, we have followed the dynamical evolution 
of star clusters by performing iV-body simulations, in order 
to catch the moment of binary formation. Upon spotting 
the first binary, we have then analyzed in detail how this 
binary formed, with how many other stars involved, and in 
what kind of step-by-step process. 

The structure of this paper is as follows. In section 2, 
we describe our method for A^-body simulations and our 
strategies for catching the moment of binary formation and 
for analyzing the event. In section 3, we show in detail how 
the binaries are formed. In section 4, we sum up. 

2. Method 

We have performed A-body simulations of star cluster 
evolution, up to and slightly past core-collapse. The method 
for the A^-body simulations is described in section 2.1. The 
timescales of binary formation are very much smaller than 
those of cluster evolution as a whole. This means that the 
need for high-resolution data around binary formation can- 
not be simply extended to a full recording of the whole clus- 
ter evolution history with the same level of detail. We have 
to be quite selective in deciding what to store and what to 
throw out. Our strategy of data compression is described 
in section 2.2. 
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2.1. N-body simulation 

We have followed the dynamical evolution of five star 
clusters with equal-mass stars without primordial binaries. 
Our first three clusters have about a thousand stars, N = 
Ik, another has TV = 4fc, and the last one has N — 16fc, 
where Ifc = 2^" = 1024. Their initial stellar distributions 
are given by Plummer's model. Stars in the three clusters 
with N = Ik are generated with different random seeds, 
which we label here as seeds 1, 2, and 3. 

For the TV-body simulations, we have used GORILLA, 
an TV-body simulation code for star clusters (Tanikawa 
& Fukushige, 2009). This code uses a fourth-order Her- 
mite scheme with individual timesteps (Makino & Aarseth, 
1992) for the time integration. The internal orbits of iso- 
lated binaries are approximated as Kepler motion. For fur- 
ther details, see TF09. 

The accuracy parameter 77, and its value at startup, 77s, 
were set to 0.01 and 0.0025, respectively. We have adopted 
the standard A^-body units, G = M = r^ = 1, (Heggie & 
Mathieu, 1986), where G is the gravitational constant, and 
M and r^ are respectively the total mass and virial radius 
of the cluster at the start of the simulation. 

We have performed the above A^-body simulations on one 
node of the FIRST cluster at the Center for Computational 
Sciences in University of Tsukuba (Umemura et al., 2008). 
Each node of the FIRST cluster is equipped with Blade- 
GRAPE, a special purpose computer for A'-body simula- 
tions (Sugimoto et al., 1990; Makino et al., 2003; Makino, 
2008). 

2.2. Strategy to catch the moment of binary formation 

We have caught the moment of binary formation manu- 
ally, for each simulation. After determining the overall evo- 
lution of each cluster by means of an A^-body simulation as 
described in the previous section, we have plotted the time 
evolution of the number of binaries, shown in figure 1. To 
be specific, we have disregarded binaries with only modest 
binding energy, since many of them dissolved again before 
getting a chance to harden indefinitely. We have chosen our 
cut-off at lOkT, where 3/2kT is the average stellar kinetic 
energy in the cluster, i.e. IkT = 1/{6N), at time i = 0, at 
the start of the simulation. The binding energy is expressed 
as 
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where 7771 and 7772 are the masses of the binary components, 
and 7'i2 and 7;i2 are the separation and relative velocity 
between them. The time in these figures is scaled by the 
initial half-mass relaxation time. The half-mass relaxation 
time is defined as 
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where rh is a half-mass radius of a cluster, and rh — 0.8 in 
the units at i = 0. 

After thus establishing the exact time of the formation 
of the first lOfcT binary, we have done a complete rerun 
of the same simulation, taking snapshots in much higher 
time resolution just before and after the moment of binary 
formation, as illustrated in figure 1. Since the binary is 
formed in or around the cluster core, it is sufficient to take 
snapshots only of the core regions. The combination of both 
temporal and spatial data reduction resulted in orders of 
magnitude savings in memory needed to store the history 
of stellar motions around the binary formation events. 

We have taken snapshots at very small intervals of 
O.Olicr.c, only 1% of the instantaneous core crossing time 



tcr.c, given by 
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where Tc is the core radius, and v^ is the stellar velocity 
dispersion in the core. For relatively large numbers of core 
stars, we have only included those in our snapshots, but 
whenever the number of stars in the core fell below a criti- 
cal value A"c,crit, we have included the Ac,crit nearest stars 
around the density center. Our choice for A"c,crit depends 
on the number of stars in the simulation: A"c.crit = 40 (A^ = 
Ifc), 100 (A^ = 4fc), and 1000 {N = 16fc). The density center 
is defined by the procedure described in Casertano & Hut 
(1985) (CH85), as modified by MHM90, in which a density 
around a given star is defined by using the distance to the 
sixth nearest star. 

For typical snapshots, most stars arc not synchronized to 
the time of that snapshot. In those cases, we have derived 
their positions and velocities through extrapolation, using 
a predictor given by Taylor series up to the first-order time 
derivatives of their accelerations, so-called jerk. 

We summarize the data of the snapshot files, one for 
each rerun, in table 1. It is clear that the file size is quite 
modest; the largest value is about 1 Gbyte, for A^ = 16fc, 
while smaller by almost two orders of magnitude for the 
other runs. Reasons for the rapid increase of file size with 
A^ are the fact that core crossing time decreases for the 
same selection of time interval in A^-body time units, and 
the increase of the number of stars that we save for larger 
A^ values. The fact that the Ifc runs have a file size that 
is not much smaller than the 4fc run stems from the more 
conservative approach with which we started the former 
runs, saving data for a longer time span. 

3. Results 

3.1. Binary formation and core evolution 

As is well known, binary formation is closely related to 
the evolution of a cluster core. When core collapse pro- 
ceeds, more than two stars interact simultaneously because 
of the high stellar density in the core. As a result of these 
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interactions, a binary is formed. In this section, we investi- 
gate the details of core evolution around the time of binary 
formation. 

Figure 2-4 show the time evolution of core radii of the 
clusters, and the moment of first binary formation. The 
core radii are derived in the same way as CH85 modified by 
MHM90, in which a density around a given star is defined 
from the sixth nearest star. In the panels with maximum 
zoom factor (zoom-2 for N = Ifc; zoom-4 for N = 4fc; and 
zoom-5 for N = 16fc) time is scaled by the instantaneous 
value of the core crossing time. We can express this 'co- 
moving time' T as 



dt 
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T effectively presents the accumulation of time on a clock 
with variable speed, in such a way that each increment in 
time is scaled by the value of the core crossing time at that 
particular moment. In this way, the very brief periods of 
deepest core collapse are guaranteed to be magnified so as 
to be completely resolved in time. The first binaries are 
formed during the time interval between the two dashed 
lines in the panels with maximum zoom factor. 

As expected, the time at which the first binaries are 
formed is comparable for the different runs, when expressed 
in A^-body time units: the first binary in each run is formed 
between ISfrh.i and 21tih,i- At the time of the formation of 
the first lOfcT binary, the core radii are 2 x 10^'^ - 2 x 10^^ 
for the A^ = Ifc runs with seed 1 and 2, 1 x 10"^ - 4 x 10"^ 
for the A^ = Ifc run with seed 3, 5 x lO'^ - 1 x 10"^ ^j. ^^^le 
N = 4k run, and 1 x 10"'' - 1 x lO'^ for the A^ = 16fc run 
(see the panels with maximum zoom factor in figure 2-4). 

When inspecting the zoom-1 panels of figures 2, 3, and 4, 
the core radii around the time pointed at by the arrows take 
on locally minimal values during intervals of O.ltrh.i before 
and after these times. This can be further investigated by 
following the time evolution of the maximum values of the 
core radii during fluctuations on smaller timescales than the 
initial half-mass relaxation time. The maximum values for 
the core radii, around the times pointed at by the arrows, 
are 2 x 10^^ for the Ifc runs, 1 x 10^^ for the 4fc run, and 
2 X 10^'^ for the 16fc run. In contrast, corresponding values 
at other times during this O.lirh.i period are 3 x 10~^ for 
the Ifc runs, 2 x 10~^ for the 4fc run, and 4 x 10"'^ for the 
16fc run. 

Looking in further detail at the core radii around the time 



of binary formation, there is no precise correlation between 
binary formation and core radii. In fact, core radii at times 
during which no binary is formed are just about as small 
as core radii at binary formation time, something that is 
especially clear in the 4fc and 16fc runs (see zoom-3 and 
zoom-4 of figure 3, and zoom-4 of figure 4, respectively). 
Core radii are almost constant during 1 x lO^^trh.i in the 
case of the 4fc run (see zoom-3 of figure 3) and 1 x 10~'^irh,i 
in the case of the 16fc run (see zoom-4 of figure 4). On the 
other hand, the binaries are formed during intervals much 
less than 2.5 x IQ-^i^j^ ^ (4fc) and 5 x lO^^ii-ha (16fc). The 
timescale on which core radii small enough to form a binary 
(about 1 X 10^^ for 4fc and 1 x 10~^ for 16fc) persist is 
several times longer than binary formation timescale. We 
are clearly dealing with events with low probability, even 
at these high densities, occurring stochastically. 

On the way toward binary formation, the core radii some- 
times become smaller than any other time by a factor of 
ten, for example r = —7 in the N = Ik run with seed 1, 
T = —7 in the N = Ik run with seed 2, r = —12 and —1 in 
the A^ = 16fc run. These decrease of core radii is not cor- 
related with binary formation. Such decrease of core radii 
occurs at the other time. 

Figure 5 shows the time evolution of the number of stars 
in the core around the time of binary formation. Note that 
the numbers of stars in the core around that time fiuctu- 
ates by roughly a factor or 10, with significant differences 
between different runs. This in itself already indicates that 
a homogeneous background density assumption for three- 
body binary formation breaks down. As we will now see in 
detail, binaries are indeed formed by processes other than 
simple three-body encounters. 



3.2. Hard Binary Formation Histories 

For more than 35 years, the conventional wisdom for 
the formation of the first hard binary in core collapse has 
postulated a relatively straightforward picture. While the 
density of the core increases, near-simultaneous three-body 
encounters become more and more likely. These can lead 
to the formation of binaries, where the third star picks up 
the excess energy, leaving the other two stars bound. Sub- 
sequent encounters are most likely to destroy these newly 
formed soft binaries, but occasionally a soft binary hardens 
against the odds, and then becomes stable when its bind- 
ing energy exceeds several fcT (Heggie, 1975; Hut, 1983; 
Goodman & Hut, 1993). 

We are now in a position to check this conventional wis- 
dom, using the information we have gathered about orbits 
of stars and the time evolution of the pair-wise distances 
between stars in the core. 

First, let us have a look at the orbits of stars involving 
binary formation in figure 6. Apparently, the binary in the 
N = Ak run, which presents an unusually simple formation 
case, is formed through an encounter among three single 
stars. The red, blue, and dark green stars simultaneously 



approach each other. Subsequently, the red and blue stars 
compose of a binary, and the dark green star is ejected from 
them. This indeed corresponds to the classical picture for 
binary formation. On the other hand, many stars involve 
the formation of the binaries in all the other runs. Their 
orbits are too complicated to determine the mechanism of 
binary formation from these figures, even though they al- 
ready suggest a mechanism quite different from the classi- 
cal picture. 

In order to find out more about what is actually hap- 
pening, we have plotted the time evolution of the pair- wise 
distances between all stars in the core, during a very small 
period of time just before binary formation. Figures 7-11 
show the interactions of stars shown in figure 6 in further 
detail. Each separate panel shows the distances of one par- 
ticular star to all other stars in the core. For each curve 
in a given panel, the color of the curve indicates the iden- 
tity of the star to which the pair-wise distance is plotted. 
When some of the curves arc intermittent, for example the 
nearest black curve in the gray panel of figure 8, it is an in- 
dication that one or another of the two stars in a pair have 
temporarily moved outside the cluster core. As described in 
section 2.2, we have not recorded the data for stars outside 
the cluster core. 

We now describe the evolution histories in detail. We will 
start with figure 7, where the results for the N = Ik run 
are given, for seed 1. Just before r = —1.5, an encounter 
involving four stars (1, 2, 4, 5) leads directly to the forma- 
tion of a hard binary (1, 5). A little later, just after r = — 1, 
that binary dissolves in a complicated dance involving five 
stars (1, 2, 3, 5, 6), and soon thereafter, around r = —0.5, 
a new hard binary (2,3) emerges from this dance. Finally, 
at r = 0, a direct three-body exchange leads to the first 
hard binary (1. 2) with an energy Ig lOkT. 

Figure 8 shows the results for the N = Ik run with seed 
2. At the start of the interval depicted, there already is a 
hard binary (8) that was formed earlier; since this binary 
is only slightly perturbed during this interval, retaining its 
identity, we have not numbered its two members separately: 
number 8 indicates both stars as a unit. Also present at the 
start of the interval is a mildly hard binary (4, 7) which 
by r = — 18 dissolves during an encounter with three other 
stars (1, 2, 5). Out of this five-body interaction emerges, 
soon thereafter around t = —17, another mildly hard bi- 
nary (1. 2). This new binary undergoes various perturba- 
tions by passing stars, hardening gradually at first. Then, 
at r = 0, a close encounter with a single star (3) leads to 
a sudden increase in binding energy, by a factor of two to 
three. 

Figure 9 shows a similar picture for the N = Ik run with 
seed 3. Around r = — 5, a three-body encounter (2, 3, 4) 
leads to the formation of a mildly hard binary (2, 3). A little 
later, at r = —3.5, a simultaneous encounter between this 
binary and two single stars (1 and 5) leads to the formation 
of a significantly harder binary (1, 2). At r = 0, a close 
encounter with a single star (3) leads to a sudden increase 
in binding energy, by about a factor of two. 



Figure 10 shows the N ~ 4fc run. In contrast to the 
preceding three cases, here hard binary formation is much 
simpler. At r = three stars (1, 2, 3) have a simultaneous 
close encounter, leading directly to the formation of the 
first hard binary (1, 2) with an energy Ig lOkT. 

Figure 11, for the N = 16fc run, is again far more com- 
phcated. At t = —17, four single stars (4, 5, 6, 8) meet, 
producing a mildly hard binary (5, 8). At r = —14 the bi- 
nary dissolves in a many-body dance involving several other 
stars that keep interacting until at r = —9a new binary (1, 
2) is formed; at the last step in its formation mainly three 
stars (1. 2. 5) were involved. Finally, at r = 0, a close en- 
counter with a single star (3) leads to a sudden increase in 
binding energy, by a factor of two to three. 

We conclude from these five cases that direct three-body 
binary formation is relatively rare. In the one case that we 
observed, in figure 10, it led to a binary with a binding 
energy of more than lOfcT, a very different picture than the 
gradual hardening of a soft binary. 

In the other four cases that we studied, the situation 
was quite a bit more complicated. In all of those cases, 
there were episodes in which at least four or five stars were 
simultaneously involved in close encounters that played a 
role in the formation of the first really hard binary. 



4. Summary 

In this paper, we have shown that core collapse of star 
clusters gives rise to the formation of hard binaries via far 
more complicated processes than have been hitherto con- 
sidered. Most importantly, the traditional approximation of 
binary formation through simultaneous encounters of three 
bodies offers a vastly oversimplified picture of what really 
happens. Typically, more stars are strongly involved, dy- 
namically, in a complex multi-body dance that in no way 
can be disentangled to lead to a hierarchy with only three 
stars in the center. 

We have found only one exception, among the five cases 
studied in great detail. Even in that case, the traditional 
picture of slow hardening of an initially quite soft binary 
doesn't apply. Rather, that three-body process of binary 
formation produced a rather hard binary very quickly, with- 
out subsequent hardening encounters with other stars. 

We can point to a single reason for the fact that the two 
core assumptions (3-body encounters, and gradual harden- 
ing) of the traditional scenario for hard binary formation 
turn out to be unrealistic. The main reason for the break- 
down of those assumptions is the very small core at the 
moment of core collapse. What we have found is that, typ- 
ically, the core at that time contains less than 10 stars, and 
often only 5 or so stars - so few that a precise measure of 
the actual number of stars becomes difficult. Given the dis- 
crete nature of core stars, the use of continuum properties 
like density for a population of only 5 stars is obviously 
rather problematic. 

If core bounce would occur at somewhat larger core sizes. 



with cores containing twenty to thirty stars or more, the 
assumptions of three-body encounters and gradual harden- 
ing may weh have been correct. Following the path break- 
ing analytical investigation of the role of binaries by Heggie 
(1975) in the seventies, subsequent analytical work in the 
eighties seems to suggest that core sizes during core bounce 
remained large enough to ensure Heggie's assumptions. For 
example, from table 2 of Goodman (1984) , equilibrium val- 
ues of core sizes between N = Ik and N = 16fc range be- 
tween twenty and sixty stars. A more detailed analysis in 
Goodman (1987) gave rise to his equation (IV. 8) suggest- 
ing a core size at deep bounce large enough for the core to 
contain dozens of stars. 

Our simulations, in contrast, suggest that deep core col- 
lapse causes the number of stars in the core to dip below 
ten in all cases, independent of the number N of stars for 
the cluster as a whole. During the final stages before the 
deepest collapse, the remainder of the cluster is effectively 
frozen, on the local time scale of the most central regions, 
and even the order of magnitude of N may become irrele- 
vant. We have demonstrated this for the range iV = Ifc to 
TV = 16k, and we predict a very similar behavior for larger 
N values, up to iV = 10'^. 

It is an interesting question, why the analytic estimates 
for the core size at the time of core bounce have overes- 
timated the number of stars in the core iVc^min by almost 
an order of magnitude. Our results indicate iVc^min ^ 5 
while the analytic estimates have predicted results more 
like A''c,min ~ 30. One possible reason could be the break- 
down of the instantaneous approximation for interactions, 
used in analytic treatments. With the timescale for forma- 
tion and subsequent hardening of binaries being compara- 
ble to the timescale of core contraction, we may need to 
extend the analytical approach. We are currently looking 
into these questions in more detail, continuing the investi- 
gations reported here. 
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Fig. 1. Time evolution of the number of binaries with more than lOfcT (nj,). Time is scaled by initial half-mass relaxation time, defined in 
equation (2). 
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Fig. 2. Time evolution of core radii in the A'^ = lA: runs with seed 1, 2, and 3. The core radii in zoom-0 panels are calculated at each 1 
A'^-body time unit (5.9 X 10~^irh). Those in the other panels are calculated at each 0.1tcr,c. The width of the time in panels zoom-1 is It^h i- 
The range of the time in panel zoom-n is the same as that between two vertical dashed lines in panel zoom-(n — 1). The ranges between two 
vertical dashed lines in panels zoom-2 are the same as those of the time in figure 7, 8, and 9. The arrows and the time r = show the time 
when the first binary is formed in each run. 
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Fig. 3. Time evolution of core radii in the N = ik run. The core radius in zoom-O panel is calculated at each 1 A'^-body time unit (1.8 X 10~'^trh). 
Those in the other panels are calculated at each O.ltcr.c The width of the time in panels zoom-1 is Itrh.i- From zoom-1 to zoom-3, the time 
is zoomed out by ten times step by step. The meanings of the dashed lines are the same as those in figure 2. The range between the two 
vertical dashed lines in panel zoom-4 is the same as that of the time in figure 10. The arrows and the time r = show the time when the 
first binary is formed. 
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Fig. 4. Time evolution of core radii in the N = 16A; run. The arrows and the time r = show the time when the first binary is formed. The 
core radius in zoom-O panel is calculated at each 1 A'^-body time unit (5.4 X 10~^tj.[j). Those in panels zoom-1, 2, 3, and 4 are calculated at 
each Itcr.c, and those in panel zoom-5 at each O.licr.c. The width of the time in panels zoom-1 is Itj-h i- From zoom-1 to zoom-4, the time 
is zoomed out by ten times step by step. The meanings of the dashed lines are the same as those in figure 2. The range between the two 
vertical dashed lines in panel zoom-4 is the same as that of the time in figure 11. 
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Fig. 5. Time evolution of the number of stars in core at around binary formation in the N = Ik runs with seed 1, 2, and 3, N = 4k run, and 
A'' = 16fe run. 
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Fig. 6. Orbits of stars involving the binary formation in the N = Ik run with seed 1. The stars whose orbits are drawn by red and blue 
curves finally become the binary components, and their directions are indicated by black arrows. The time ranges are the same as those in 
panels zoom-2, 4, and 5 in figure 2, 3, and 4, respectively. 
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Fig. 7. Time evolution of separations among stars involving binary formation in the A^ = lA: run with seed 1. The colors are the same as 
those in figure 6. The black curves indicate stars which do not involve the binary formation. For each star, following the color, the number 
between parentheses refers to the numbers used in the discussion in section 3.2. 
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Fig. 8. Time evolution of separations among stars involving binary formation in the N = Ik run with seed 2. 
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Fig. 9. Time evolution of separations among stars involving binary formation in the N = Ik run with seed 3. 



o 



o 



o 




Fig. 10. Time evolution of separations among stars involving binary formation in the N = ik run. 
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Fig. 11. Time evolution of separations among stars involving binary formation in the A'^ = 16fc run. 



